

d=read_ascii('strat.dat')
r=reform(d.field1[0,*])
t=d.field1[1,*]
rho=d.field1[2,*]
rho_ad=d.field1[3,*]
temp=d.field1[4,*]
th=d.field1[5,*]
p=rho*t
dth=deriv(r,th)

ll=(size(rho))[2]
print,ll
theta=!pi*indgen(ll)/(ll-1)


GG=6.674e-11            ;N(m/kg)^2
solar_mass=1.988e30     ;kg
solar_radii=6.96e8
gamma=5./3

;g=GG*solar_mass/(solar_radii*r)^2
;omega=2.*!pi/(28.*24.*3600.)

;N2 = g/th*dth
;BV=r*(deriv(r,p)/(gamma*p)-deriv(r,rho)/rho)
yr=[min(th),max(th)]
plot,r,th,yr=yr,ystyle=1
;oplot,[0,1],[th00,th00],li=2

print,'rho_max/rho_min =',max(rho)/min(rho)
print,"\Delta\Theta1",max(th)-th[0]
print,"\Delta\Theta2",max(th)-th[ll-1]
jj=where(r ge 0.73)
thup=th[jj]
rup=r[jj]
yrup=[min(thup),max(thup)]

;th2d=fltarr(ll,ll)
;
;for i=0,ll-1 do begin
;  for j=0, ll-1 do begin
;;    th2d[i,j]=th[i]+0.5*(3.*cos(2.*theta[j])-1)
;    th2d[i,j]=1.*Legendre(cos(theta[j]),2);-.1*Legendre(cos(theta[j]),4)
;  endfor
;endfor

;x=r#sin(theta)
;y=r#cos(theta)

;contour,th2d,x,y,/fi,nl=64,/iso
;print,"lat_dif= ",max(th2d[127,*])-min(th2d[127,*])


end
